#pragma once

using namespace System;
using namespace System::IO;
using namespace NUnit::Framework;
using namespace LatoolNet;

namespace LatoolNetTest {

	[TestFixture]
	public ref class EigenproblemTest {
	public:

		[Test]
		void TestSymmetricEigenproblem() {

			int rownum = 4;
			int colnum = 4;
			Matrix^ a = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);

			a[0, 0] = 1;
			a[0, 1] = 2;
			a[0, 2] = 3;
			a[0, 3] = 4;

			a[1, 1] = 2;
			a[1, 2] = 3;
			a[1, 3] = 4;

			a[2, 2] = 3;
			a[2, 3] = 4;

			a[3, 3] = 4;

			Matrix ^ orig = a->Clone();

			array<double>^ resultValues;
			Matrix^ resultVectors;

			Eigenproblem::Solve(a, resultValues, resultVectors);

			for (int n = 0; n < resultValues->Length; n++) {
				double lambda = resultValues[n];

				Matrix^ vector = gcnew Matrix(4, 1);
				for (int i = 0; i < 4; i++) {
					vector[i, 0] = resultVectors[i, n];
				}

				Matrix^ lhs = orig * vector;
				Matrix^ rhs = lambda * vector;

				for (int i = 0; i < 4; i++) {
					Assert::AreEqual(lhs[i, 0], rhs[i, 0], 1e-10, "Test: EigenvalueProblem");
				}
			}
		};

		[Test]
		void TestSymmetricEigenproblem2() {

			int rownum = 4;
			int colnum = 4;
			Matrix^ a = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);

			a[0, 0] = 1;
			a[0, 1] = 2;
			a[0, 2] = 3;
			a[0, 3] = 4;

			a[1, 1] = 2;
			a[1, 2] = 3;
			a[1, 3] = 4;

			a[2, 2] = 3;
			a[2, 3] = 4;

			a[3, 3] = 4;

			array<double>^ resultValues1;
			Matrix^ resultVectors1;

			Eigenproblem::Solve(a->Clone(), resultValues1, resultVectors1);

			array<double>^ resultValues2;
			Eigenproblem::Solve(a->Clone(), resultValues2);

			for (int i = 0; i < 4; i++) {
				Assert::AreEqual(resultValues1[i], resultValues2[i], 1e-10, "Test: EigenvalueProblem");
			}
		};


		[Test]
		void TestHermitianEigenproblem() {

			int size = 5;
			int rownum = size;
			int colnum = size;
			ComplexMatrix^ a = gcnew ComplexMatrix(rownum, colnum, MatrixType::ComplexHermitian);

			a[0, 0] = Complex(1.0, 0.0);
			a[0, 1] = Complex(2.0, -1.0);
			a[0, 2] = Complex(3.0, -1.0);
			a[0, 3] = Complex(0.0, 0.0);
			a[0, 4] = Complex(0.0, 0.0);

			a[1, 1] = Complex(2.0, 0.0);
			a[1, 2] = Complex(3.0, -2.0);
			a[1, 3] = Complex(4.0, -2.0);
			a[1, 4] = Complex(0.0, 0.0);

			a[2, 2] = Complex(3.0, 0.0);
			a[2, 3] = Complex(4.0, -3.0);
			a[2, 4] = Complex(5.0, -3.0);

			a[3, 3] = Complex(4.0, 0.0);
			a[3, 4] = Complex(5.0, -4.0);

			a[4, 4] = Complex(5.0, 0.0);

			ComplexMatrix ^ orig = a->Clone();

			array<double>^ resultValues;
			ComplexMatrix^ resultVectors;

			Eigenproblem::Solve(a, resultValues, resultVectors);

			for (int n = 0; n < resultValues->Length; n++) {
				double lambda = resultValues[n];

				ComplexMatrix^ vector = gcnew ComplexMatrix(size, 1);
				for (int i = 0; i < size; i++) {
					vector[i, 0] = resultVectors[i, n];
				}

				ComplexMatrix^ lhs = orig * vector;
				ComplexMatrix^ rhs = lambda * vector;

				for (int i = 0; i < 4; i++) {
					Assert::AreEqual(lhs[i, 0].Real, rhs[i, 0].Real, 1e-10, "Test: EigenvalueProblem");
					Assert::AreEqual(lhs[i, 0].Imag, rhs[i, 0].Imag, 1e-10, "Test: EigenvalueProblem");
				}
			}
		};

		[Test]
		void TestHermitianEigenproblem2() {

			int size = 5;
			int rownum = size;
			int colnum = size;
			ComplexMatrix^ a = gcnew ComplexMatrix(rownum, colnum, MatrixType::ComplexHermitian);

			a[0, 0] = Complex(1.0, 0.0);
			a[0, 1] = Complex(2.0, -1.0);
			a[0, 2] = Complex(3.0, -1.0);
			a[0, 3] = Complex(0.0, 0.0);
			a[0, 4] = Complex(0.0, 0.0);

			a[1, 1] = Complex(2.0, 0.0);
			a[1, 2] = Complex(3.0, -2.0);
			a[1, 3] = Complex(4.0, -2.0);
			a[1, 4] = Complex(0.0, 0.0);

			a[2, 2] = Complex(3.0, 0.0);
			a[2, 3] = Complex(4.0, -3.0);
			a[2, 4] = Complex(5.0, -3.0);

			a[3, 3] = Complex(4.0, 0.0);
			a[3, 4] = Complex(5.0, -4.0);

			a[4, 4] = Complex(5.0, 0.0);

			array<double>^ resultValues1;
			ComplexMatrix^ resultVectors1;

			Eigenproblem::Solve(a->Clone(), resultValues1, resultVectors1);

			array<double>^ resultValues2;
			Eigenproblem::Solve(a->Clone(), resultValues2);

			for (int i = 0; i < 5; i++) {
				Assert::AreEqual(resultValues1[i], resultValues2[i], 1e-10, "Test: EigenvalueProblem");
			}
		};

		[Test]
		void TestGenralizedSymmetricDefiniteEigenproblem() {

			int rownum = 4;
			int colnum = 4;
			Matrix^ a = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);

			a[0, 0] = 0.24;
			a[0, 1] = 0.39;
			a[0, 2] = 0.42;
			a[0, 3] = -0.16;

			a[1, 1] = -0.11;
			a[1, 2] = 0.79;
			a[1, 3] = 0.63;

			a[2, 2] = -0.25;
			a[2, 3] = 0.48;

			a[3, 3] = -0.03;

			Matrix^ b = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);
			b[0, 0] = 4.16;
			b[0, 1] = -3.12;
			b[0, 2] = 0.56;
			b[0, 3] = -0.10;

			b[1, 1] = 5.03;
			b[1, 2] = -0.83;
			b[1, 3] = 1.09;

			b[2, 2] = 0.76;
			b[2, 3] = 0.34;

			b[3, 3] = 1.18;

			Matrix ^ origa = a->Clone();
			Matrix ^ origb = b->Clone();

			array<double>^ resultValues;
			Matrix^ resultVectors;

			Eigenproblem::Solve(a, b, resultValues, resultVectors);

			for (int n = 0; n < resultValues->Length; n++) {
				double lambda = resultValues[n];
				Matrix^ vector = gcnew Matrix(4, 1);
				for (int i = 0; i < 4; i++) {
					vector[i, 0] = resultVectors[i, n];
				}

				Matrix^ lhs = origa * vector;
				Matrix^ rhs = lambda * origb * vector;

				for (int i = 0; i < 4; i++) {
					Assert::AreEqual(lhs[i, 0], rhs[i, 0], 1e-10, "Test: EigenvalueProblem");
				}
			}
		};

		[Test]
		void TestGenralizedSymmetricDefiniteEigenproblem2() {

			int rownum = 4;
			int colnum = 4;
			Matrix^ a = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);

			a[0, 0] = 0.24;
			a[0, 1] = 0.39;
			a[0, 2] = 0.42;
			a[0, 3] = -0.16;

			a[1, 1] = -0.11;
			a[1, 2] = 0.79;
			a[1, 3] = 0.63;

			a[2, 2] = -0.25;
			a[2, 3] = 0.48;

			a[3, 3] = -0.03;

			Matrix^ b = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);
			b[0, 0] = 4.16;
			b[0, 1] = -3.12;
			b[0, 2] = 0.56;
			b[0, 3] = -0.10;

			b[1, 1] = 5.03;
			b[1, 2] = -0.83;
			b[1, 3] = 1.09;

			b[2, 2] = 0.76;
			b[2, 3] = 0.34;

			b[3, 3] = 1.18;

			array<double>^ resultValues1;
			Matrix^ resultVectors1;

			Eigenproblem::Solve(a->Clone(), b->Clone(), resultValues1, resultVectors1);

			array<double>^ resultValues2;

			Eigenproblem::Solve(a->Clone(), b->Clone(), resultValues2);

			for (int i = 0; i < 4; i++) {
				Assert::AreEqual(resultValues1[i], resultValues2[i], 1e-10, "Test: EigenvalueProblem");
			}
		};

		[Test]
		void TestGenralizedHermitianDefiniteEigenproblem() {

			int rownum = 4;
			int colnum = 4;
			ComplexMatrix^ a = gcnew ComplexMatrix(rownum, colnum, MatrixType::ComplexHermitian);

			a[0, 0] = Complex(-7.36, 0.0);
			a[0, 1] = Complex(0.77, -0.43);
			a[0, 2] = Complex(-0.64, -0.92);
			a[0, 3] = Complex(3.01, -6.97);

			a[1, 1] = Complex(3.49, 0.0);
			a[1, 2] = Complex(2.19, 4.45);
			a[1, 3] = Complex(1.90, 3.73);

			a[2, 2] = Complex(0.12, 0.0);
			a[2, 3] = Complex(2.88, -3.17);

			a[3, 3] = Complex(-2.54, 0.0);

			ComplexMatrix ^ origa = a->Clone();

			ComplexMatrix^ b = gcnew ComplexMatrix(rownum, colnum, MatrixType::ComplexHermitian);

			b[0, 0] = Complex(3.23, 0.0);
			b[0, 1] = Complex(1.51, -1.92);
			b[0, 2] = Complex(1.90, 0.84);
			b[0, 3] = Complex(0.42, 2.50);

			b[1, 1] = Complex(3.58, 0.0);
			b[1, 2] = Complex(-0.23, 1.11);
			b[1, 3] = Complex(-1.18, 1.37);

			b[2, 2] = Complex(4.09, 0.0);
			b[2, 3] = Complex(2.33, -0.14);

			b[3, 3] = Complex(4.29, 0.0);

			ComplexMatrix^ origb = b->Clone();

			array<double>^ resultValues;
			ComplexMatrix^ resultVectors;

			Eigenproblem::Solve(a, b, resultValues, resultVectors);

			for (int n = 0; n < resultValues->Length; n++) {
				double lambda = resultValues[n];
				ComplexMatrix^ vector = gcnew ComplexMatrix(4, 1);
				for (int i = 0; i < 4; i++) {
					vector[i, 0] = resultVectors[i, n];
				}

				ComplexMatrix^ lhs = origa * vector;
				ComplexMatrix^ rhs = lambda * origb * vector;

				for (int i = 0; i < 4; i++) {
					Assert::AreEqual(lhs[i, 0].Real, rhs[i, 0].Real, 1e-10, "Test: EigenvalueProblem");
					Assert::AreEqual(lhs[i, 0].Imag, rhs[i, 0].Imag, 1e-10, "Test: EigenvalueProblem");

				}
			}
		};

		[Test]
		void TestGenralizedHermitianDefiniteEigenproblem2() {

			int rownum = 4;
			int colnum = 4;
			ComplexMatrix^ a = gcnew ComplexMatrix(rownum, colnum, MatrixType::ComplexHermitian);

			a[0, 0] = Complex(-7.36, 0.0);
			a[0, 1] = Complex(0.77, -0.43);
			a[0, 2] = Complex(-0.64, -0.92);
			a[0, 3] = Complex(3.01, -6.97);

			a[1, 1] = Complex(3.49, 0.0);
			a[1, 2] = Complex(2.19, 4.45);
			a[1, 3] = Complex(1.90, 3.73);

			a[2, 2] = Complex(0.12, 0.0);
			a[2, 3] = Complex(2.88, -3.17);

			a[3, 3] = Complex(-2.54, 0.0);

			ComplexMatrix^ b = gcnew ComplexMatrix(rownum, colnum, MatrixType::ComplexHermitian);

			b[0, 0] = Complex(3.23, 0.0);
			b[0, 1] = Complex(1.51, -1.92);
			b[0, 2] = Complex(1.90, 0.84);
			b[0, 3] = Complex(0.42, 2.50);

			b[1, 1] = Complex(3.58, 0.0);
			b[1, 2] = Complex(-0.23, 1.11);
			b[1, 3] = Complex(-1.18, 1.37);

			b[2, 2] = Complex(4.09, 0.0);
			b[2, 3] = Complex(2.33, -0.14);

			b[3, 3] = Complex(4.29, 0.0);

			array<double>^ resultValues1;
			ComplexMatrix^ resultVectors1;

			Eigenproblem::Solve(a->Clone(), b->Clone(), resultValues1, resultVectors1);

			array<double>^ resultValues2;

			Eigenproblem::Solve(a->Clone(), b->Clone(), resultValues2);

			for (int i = 0; i < 4; i++) {
				Assert::AreEqual(resultValues1[i], resultValues2[i], 1e-10, "Test: EigenvalueProblem");

			}
		};

	};
}
